### RLM MODEL RESULTS ###

import numpy as np

def _shift_intercept(arr):
    """
    A convenience function to make the SAS covariance matrix
    compatible with statsmodels.rlm covariance
    """
    arr = np.asarray(arr)
    side = int(np.sqrt(len(arr)))
    return np.roll(np.roll(arr.reshape(side,side),-1, axis =1), -1, axis=0)

class Huber(object):
    """
    """
    def __init__(self):
        self.params = np.array([  0.82937387,   0.92610818,  -0.12784916,
            -41.02653105])
        self.bse = np.array([ 0.11118035,  0.3034081 ,  0.12885259,  9.8073472 ])
        self.scale = 2.4407137948148447
        self.weights = np.array([ 1.,  1.,  0.7858871,  0.50494094,  1.,
            1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
            0.36814106])
        self.resid  = np.array([ 3.05027584, -2.07757332,  4.17721071,
            6.50163171, -1.64615192, -2.57226011, -1.73127333, -0.73127333,
            -2.25476463,  0.48083217,1.63147461,  1.42973363, -2.26346951,
            -0.78323693,  2.26646556, 0.88291808, -0.83307835,  0.06186577,
            0.26360675,  1.54306186, -8.91752986])
        self.df_model = 3
        self.df_resid = 17
        self.bcov_unscaled = np.array([[  1.72887367e-03,  -3.47079127e-03,
             -6.79080082e-04, 2.73387119e-02],
           [ -3.47079127e-03,   1.28754242e-02,   9.95952051e-07,
             -6.19611175e-02],
           [ -6.79080082e-04,   9.95952051e-07,   2.32216722e-03,
             -1.59355028e-01],
           [  2.73387119e-02,  -6.19611175e-02,  -1.59355028e-01,
              1.34527267e+01]]) # From R
        self.fittedvalues = np.array([ 38.94972416,  39.07757332,  32.82278929,
            21.49836829,
            19.64615192,  20.57226011,  20.73127333,  20.73127333,
            17.25476463,  13.51916783,  12.36852539,  11.57026637,
            13.26346951,  12.78323693,   5.73353444,   6.11708192,
             8.83307835,   7.93813423,   8.73639325,  13.45693814,
             23.91752986])
        self.tvalues = np.array([ 7.45971657,  3.0523516 , -0.99221261,
            -4.18324448])
        # from R this is equivalent to
#        self.res1.params/np.sqrt(np.diag(self.res1.bcov_scaled))

#    def conf_int(self): # method to be consistent with sm
#        return



    # The below are taken from SAS

    huber_h1 = [95.8813, 0.19485, -0.44161, -1.13577, 0.1949, 0.01232,
            -0.02474, -0.00484, -0.4416, -0.02474, 0.09177, 0.00001, -1.1358,
            -0.00484, 0.00001, 0.01655]
    h1 = _shift_intercept(huber_h1)

    huber_h2 = [82.6191, 0.07942, -0.23915, -0.95604, 0.0794, 0.01427,
            -0.03013, -0.00344, -0.2392, -0.03013, 0.10391, -0.00166, -0.9560,
            -0.00344, -0.00166, 0.01392]
    h2 = _shift_intercept(huber_h2)

    huber_h3 = [70.1633, -0.04533, -0.00790, -0.78618, -0.0453, 0.01656,
            -0.03608, -0.00203, -0.0079, -0.03608,  0.11610, -0.00333, -0.7862,
            -0.00203, -0.00333,  0.01138]
    h3 = _shift_intercept(huber_h3)

class Hampel(object):
    """
    """
    def __init__(self):
        self.params = np.array([  0.74108304,   1.22507934,  -0.14552506,
            -40.47473236])
        self.bse = np.array([  0.13482596,   0.36793632,   0.1562567 ,
            11.89315426])
        self.scale = 3.0882646556217064
        self.weights = np.array([ 1.,  1.,  1.,  1.,  1.,   1.,  1.,  1.,  1.,
            1., 1.,  1.,  1.,  1.,  1., 1.,  1.,  1.,  1.,  1., 0.80629719])
        self.resid = np.array([ 3.06267708, -2.08284798,  4.36377602,
            5.78635972, -1.7634816 ,
           -2.98856094, -2.34048993, -1.34048993, -3.02422878,  1.08249252,
            2.39221804,  2.47177232, -1.62645737, -0.25076107,  2.32088237,
            0.88430719, -1.37812296, -0.35944755, -0.43900184,  1.40555003,
           -7.65988702])
        self.df_model = 3
        self.df_resid = 17
        self.bcov_unscaled = np.array([[  1.72887367e-03,  -3.47079127e-03,
            -6.79080082e-04,
              2.73387119e-02],
           [ -3.47079127e-03,   1.28754242e-02,   9.95952051e-07,
             -6.19611175e-02],
           [ -6.79080082e-04,   9.95952051e-07,   2.32216722e-03,
             -1.59355028e-01],
           [  2.73387119e-02,  -6.19611175e-02,  -1.59355028e-01,
              1.34527267e+01]])
        self.fittedvalues = np.array([ 38.93732292,  39.08284798,
            32.63622398,  22.21364028,
            19.7634816 ,  20.98856094,  21.34048993,  21.34048993,
            18.02422878,  12.91750748,  11.60778196,  10.52822768,
            12.62645737,  12.25076107,   5.67911763,   6.11569281,
             9.37812296,   8.35944755,   9.43900184,  13.59444997,
             22.65988702])
        self.tvalues = np.array([ 5.49659011,  3.32959607, -0.93132046,
            -3.40319578])

    hampel_h1 = [141.309,  0.28717, -0.65085, -1.67388, 0.287,  0.01816,
            -0.03646, -0.00713, -0.651, -0.03646,  0.13524,  0.00001, -1.674,
            -0.00713, 0.00001,  0.02439]
    h1 = _shift_intercept(hampel_h1)

    hampel_h2 = [135.248,  0.18207, -0.36884, -1.60217, 0.182, 0.02120,
            -0.04563, -0.00567, -0.369, -0.04563,  0.15860, -0.00290, -1.602,
            -0.00567, -0.00290, 0.02329]
    h2 = _shift_intercept(hampel_h2)

    hampel_h3 = [128.921,  0.05409, -0.02445, -1.52732, 0.054,  0.02514,
            -0.05732, -0.00392, -0.024, -0.05732,  0.18871, -0.00652, -1.527,
            -0.00392, -0.00652,  0.02212]
    h3 = _shift_intercept(hampel_h3)

class BiSquare(object):
    def __init__(self):
        self.params = np.array([  0.9275471 ,   0.65073222,  -0.11233103,
            -42.28525369])
        self.bse = np.array([ 0.10805398,  0.29487634,  0.12522928,
            9.5315672 ])
        self.scale = 2.2818858795649497
        self.weights = np.array([ 0.89283149,  0.88496132,  0.79040651,
            0.3358111 ,  0.94617358, 0.90040725,  0.96630596,  0.99729171,
            0.94968061,  0.99900087, 0.98959903,  0.99831448,  0.84731833,
            0.96455873,  0.91767906,  0.98724523,  0.99762848,  0.99694419,
            0.98650731,  0.95897484, 0.00222999])
        self.resid = np.array([  2.50917802,  -2.60315301,   3.56070896,
            6.93256033, -1.76597524,  -2.41670746,  -1.39345348,  -0.39345348,
            -1.70651907,  -0.23917521,   0.77180408,   0.31020526,
            -3.01451315,  -1.42960401,   2.19218084,   0.85518774,
            -0.36817892,   0.4181383 ,   0.87973712,   1.53911661,
            -10.43556344])
        self.df_model = 3
        self.df_resid = 17
        self.bcov_unscaled = np.array([[  1.72887367e-03,  -3.47079127e-03,
             -6.79080082e-04, 2.73387119e-02],
            [ -3.47079127e-03,   1.28754242e-02,   9.95952051e-07,
             -6.19611175e-02],
           [ -6.79080082e-04,   9.95952051e-07,   2.32216722e-03,
             -1.59355028e-01],
           [  2.73387119e-02,  -6.19611175e-02,  -1.59355028e-01,
              1.34527267e+01]])
        self.fittedvalues = np.array([ 39.49082198,  39.60315301,  33.43929104,
            21.06743967,
            19.76597524,  20.41670746,  20.39345348,  20.39345348,
            16.70651907,  14.23917521,  13.22819592,  12.68979474,
            14.01451315,  13.42960401,   5.80781916,   6.14481226,
             8.36817892,   7.5818617 ,   8.12026288,  13.46088339,
             25.43556344])
        self.tvalues = np.array([ 8.58410823,  2.20679698, -0.8970029 ,
            -4.43633799])


    bisquare_h1 = [90.3354,  0.18358, -0.41607, -1.07007, 0.1836, 0.01161,
            -0.02331, -0.00456, -0.4161, -0.02331,  0.08646, 0.00001, -1.0701,
            -0.00456, 0.00001,  0.01559]
    h1 = _shift_intercept(bisquare_h1)

    bisquare_h2 = [67.82521, 0.091288, -0.29038, -0.78124, 0.091288,
            0.013849, -0.02914, -0.00352, -0.29038, -0.02914, 0.101088, -0.001,
            -0.78124, -0.00352,   -0.001, 0.011766]
    h2 = _shift_intercept(bisquare_h2)

    bisquare_h3 = [48.8983, 0.000442, -0.15919, -0.53523, 0.000442,
            0.016113, -0.03461, -0.00259, -0.15919, -0.03461, 0.112728,
            -0.00164, -0.53523, -0.00259, -0.00164, 0.008414]
    h3 = _shift_intercept(bisquare_h3)

class Andrews(object):
    def __init__(self):
        self.params = [0.9282, 0.6492, -.1123,-42.2930]
        self.bse = [.1061, .2894, .1229, 9.3561]
        self.scale = 2.2801
        self.df_model = 3.
        self.df_resid = 17.
#        self.bcov_unscaled = [] # not given as part of SAS
        self.resid = [2.503338458, -2.608934536, 3.5548678338, 6.9333705014,
            -1.768179527, -2.417404513, -1.392991531, -0.392991531,
            -1.704759385,-0.244545418, 0.7659115325, 0.3028635237,
            -3.019999429,-1.434221475,2.1912017882, 0.8543828047,
            -0.366664104,0.4192468573,0.8822948661,1.5378731634,
            -10.44592783]

        self.sresids = [1.0979293816, -1.144242351, 1.5591155202, 3.040879735,
            -0.775498914, -1.06023995, -0.610946684, -0.172360612,
            -0.747683723, -0.107254214, 0.3359181307, 0.1328317233,
            -1.324529688, -0.629029563, 0.9610305856, 0.3747203984,
            -0.160813769, 0.1838758324, 0.3869622398, 0.6744897502,
            -4.581438458]

        self.weights = [0.8916509101, 0.8826581922, 0.7888664106, 0.3367252734,
            0.9450252405, 0.8987321912, 0.9656622, 0.9972406688,
            0.948837669, 0.9989310017, 0.9895434667, 0.998360628,
            0.8447116551, 0.9636222149, 0.916330067, 0.9869982597,
            0.9975977354, 0.9968600162, 0.9861384742, 0.9582432444, 0]

    def conf_int(self): # method to be consistent with sm
        return [(0.7203,1.1360),(.0819,1.2165),(-.3532,.1287),
                (-60.6305,-23.9555)]



    andrews_h1 = [87.5357, 0.177891, -0.40318, -1.03691, 0.177891,  0.01125,
            -0.02258, -0.00442, -0.40318, -0.02258, 0.083779, 6.481E-6,
            -1.03691, -0.00442, 6.481E-6,  0.01511]
    h1 = _shift_intercept(andrews_h1)

    andrews_h2 = [66.50472,  0.10489,  -0.3246, -0.76664, 0.10489, 0.012786,
            -0.02651,  -0.0036, -0.3246, -0.02651,  0.09406, -0.00065,
            -0.76664,  -0.0036, -0.00065, 0.011567]
    h2 = _shift_intercept(andrews_h2)

    andrews_h3 = [48.62157, 0.034949, -0.24633, -0.53394, 0.034949, 0.014088,
                -0.02956, -0.00287, -0.24633, -0.02956, 0.100628, -0.00104,
                -0.53394, -0.00287, -0.00104, 0.008441]
    h3 = _shift_intercept(andrews_h3)




### RLM Results with Huber's Proposal 2 ###
### Obtained from SAS ###

class HuberHuber(object):
    def __init__(self):
        self.h1 = [114.4936, 0.232675, -0.52734, -1.35624, 0.232675, 0.014714,
                -0.02954, -0.00578, -0.52734, -0.02954, 0.10958, 8.476E-6,
                -1.35624, -0.00578, 8.476E-6, 0.019764]
        self.h1 = _shift_intercept(self.h1)
        self.h2 = [103.2876, 0.152602, -0.33476, -1.22084, 0.152602, 0.016904,
                -0.03766, -0.00434, -0.33476, -0.03766, 0.132043, -0.00214,
                -1.22084, -0.00434, -0.00214, 0.017739]
        self.h2 = _shift_intercept(self.h2)
        self.h3 = [ 91.7544, 0.064027, -0.11379, -1.08249, 0.064027, 0.019509,
                -0.04702, -0.00278, -0.11379, -0.04702, 0.157872, -0.00462,
                -1.08249, -0.00278, -0.00462, 0.015677]
        self.h3 = _shift_intercept(self.h3)
        self.resid = [2.909155172, -2.225912162, 4.134132661, 6.163172632,
                -1.741815737, -2.789321552, -2.02642336, -1.02642336,
                -2.593402734, 0.698655, 1.914261011, 1.826699492, -2.031210331,
                -0.592975466, 2.306098648, 0.900896645, -1.037551854,
                -0.092080512, -0.004518993, 1.471737448, -8.498372406]
        self.sresids = [0.883018497, -0.675633129, 1.25483702, 1.870713355,
                -0.528694904, -0.84664529, -0.615082113, -0.311551209,
                -0.787177874, 0.212063383, 0.581037374, 0.554459746,
                -0.616535106, -0.179986379, 0.699972205, 0.273449972,
                -0.314929051, -0.027949281, -0.001371654, 0.446717797,
                -2.579518651]
        self.weights = [1, 1, 1, 0.718977066, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 1, 1, 1, 0.52141511]

        self.params = (.7990,1.0475,-0.1351,-41.0892)
        self.bse = (.1213,.3310,.1406,10.7002)
        self.scale = 3.2946
        self.df_model = 3
        self.df_resid = 17

    def conf_int(self): # method for consistency with sm
        return [(0.5612,1.0367),(.3987,1.6963),
                    (-.4106,.1405),(-62.0611,-20.1172)]



class HampelHuber(object):
    def __init__(self):
        self.h1 = [147.4727, 0.299695, -0.67924, -1.7469, 0.299695, 0.018952,
                -0.03805, -0.00744, -0.67924, -0.03805, 0.141144, 0.000011,
                -1.7469, -0.00744, 0.000011, 0.025456]
        self.h1 = _shift_intercept(self.h1)
        self.h2 = [141.148, 0.190007, -0.38493, -1.67206, 0.190007, 0.02213,
                -0.04762, -0.00592, -0.38493, -0.04762, 0.165518, -0.00303,
                -1.67206, -0.00592, -0.00303, 0.024301]
        self.h2 = _shift_intercept(self.h2)
        self.h3 = [134.5444, 0.05645, -0.02552, -1.59394, 0.05645, 0.026232,
                -0.05982, -0.00409, -0.02552, -0.05982, 0.196946, -0.0068,
                -1.59394, -0.00409, -0.0068, 0.023083]
        self.h3 = _shift_intercept(self.h3)
        self.resid = [3.125725599, -2.022218392, 4.434082972, 5.753880172,
                -1.744479058, -2.995299443, -2.358455878, -1.358455878,
                -3.068281354, 1.150212629, 2.481708553, 2.584584946,
                -1.553899388, -0.177335865, 2.335744732, 0.891912757,
                -1.43012351, -0.394515569, -0.497391962, 1.407968887,
                -7.505098501]
        self.sresids = [0.952186413, -0.616026205, 1.350749906, 1.752798302,
                -0.531418771, -0.912454834, -0.718453867, -0.413824947,
                -0.934687235, 0.350388031, 0.756000196, 0.787339321,
                -0.473362692, -0.054021633, 0.711535395, 0.27170242,
                -0.43565698, -0.120180852, -0.151519976, 0.428908041,
                -2.28627005]
        self.weights = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 1, 0.874787298]

        self.params = (.7318,1.2508,-0.1479,-40.2712)
        self.bse = (.1377, .3757, .1596, 12.1438)
        self.scale = 3.2827
        self.df_model = 3
        self.df_resid = 17

    def conf_int(self):
        return [(0.4619,1.0016),(.5145,1.9872),
                (-.4607,.1648),(-64.0727,-16.4697)]

class BisquareHuber(object):
    def __init__(self):
        self.h1 = [129.9556, 0.264097, -0.59855, -1.5394, 0.264097,
                0.016701, -0.03353, -0.00656, -0.59855, -0.03353,
                0.124379, 9.621E-6, -1.5394, -0.00656, 9.621E-6, 0.022433]
        self.h1 = _shift_intercept(self.h1)
        self.h2 = [109.7685, 0.103038, -0.25926, -1.28355, 0.103038, 0.0214,
                -0.04688, -0.00453, -0.25926, -0.04688, 0.158535, -0.00327,
                -1.28355, -0.00453, -0.00327, 0.018892]
        self.h2 = _shift_intercept(self.h2)
        self.h3 = [91.80527, -0.09171, 0.171716, -1.05244, -0.09171,
                0.027999, -0.06493, -0.00223, 0.171716, -0.06493, 0.203254,
                -0.0071, -1.05244, -0.00223, -0.0071, 0.015584]
        self.h3 = _shift_intercept(self.h3)
        self.resid = [3.034895447, -2.09863887, 4.229870063, 6.18871385,
            -1.715906134, -2.763596142, -2.010080245, -1.010080245,
            -2.590747917, 0.712961901, 1.914770759, 1.82892645, -2.019969464,
            -0.598781979, 2.260467209, 0.859864256, -1.057306197, -0.122565974,
            -0.036721665, 1.471074632, -8.432085298]
        self.sresids = [0.918227061, -0.634956635, 1.279774287, 1.872435025,
            -0.519158394, -0.836143718, -0.608162656, -0.305606249, -0.78384738,                     0.215711191, 0.579326161, 0.553353415, -0.611154703, -0.181165324,
            0.683918836, 0.26015744, -0.319894764, -0.037083121, -0.011110375,
            0.445083055, -2.551181429]
        self.weights = [0.924649089, 0.963600796, 0.856330585, 0.706048833,
            0.975591792, 0.937309703, 0.966582366, 0.991507994, 0.944798311,
            0.995764589, 0.969652425, 0.972293856, 0.966255569, 0.997011618,
            0.957833493, 0.993842376, 0.990697247, 0.9998747, 0.999988752,
            0.982030803, 0.494874977]
        self.params = (.7932, 1.0477, -0.1335, -40.8949)
        self.bse = (.1292, .3527, .1498, 11.3998)
        self.scale = 3.3052
        self.df_model = 3
        self.df_resid = 17

    def conf_int(self):
        return [(0.5399,1.0465),(.3565,1.7389),
                (-.4271,.1600),(-63.2381,-18.5517)]



class AndrewsHuber(object):
    def __init__(self):
        self.h1 = [129.9124, 0.264009, -0.59836, -1.53888, 0.264009,
                0.016696, -0.03352, -0.00656, -0.59836, -0.03352, 0.124337,
                9.618E-6, -1.53888, -0.00656, 9.618E-6, 0.022425]
        self.h1 = _shift_intercept(self.h1)
        self.h2 = [109.7595, 0.105022, -0.26535, -1.28332, .105022, 0.021321,
                -0.04664, -0.00456, -0.26535, -0.04664, 0.157885, -0.00321,
                -1.28332, -0.00456, -0.00321, 0.018895]
        self.h2 = _shift_intercept(self.h2)
        self.h3 = [91.82518, -0.08649, 0.155965, -1.05238, -0.08649, 0.027785,
                -0.06427, -0.0023, 0.155965, -0.06427, 0.201544, -0.00693,
                -1.05238, -0.0023, -0.00693, 0.015596]
        self.h3 = _shift_intercept(self.h3)
        self.resid = [3.040515104, -2.093093543, 4.235081748, 6.188729166,
                -1.714119676, -2.762695255, -2.009618953, -1.009618953,
                -2.591649784, 0.715967584, 1.918445405, 1.833412337,
                -2.016815123, -0.595695587, 2.260536347, 0.859710406,
                -1.059386228, -0.1241257, -0.039092633, 1.471556455,
                -8.424624872]
        self.sresids = [0.919639919, -0.633081011, 1.280950793, 1.871854667,
                -0.518455862, -0.835610004, -0.607833129, -0.305371248,
                -0.783875269, 0.216552902, 0.580256606, 0.554537345,
                -0.610009696, -0.180175208, 0.683726076, 0.260029627,
                -0.320423952, -0.037543293, -0.011824031, 0.445089734,
                -2.548127888]
        self.weights = [0.923215335, 0.963157359, 0.854300342, 0.704674258,
                0.975199805, 0.936344742, 0.9660077, 0.991354016, 0.943851708,
                0.995646409, 0.968993767, 0.971658421, 0.965766352, 0.99698502,
                0.957106815, 0.993726436, 0.990483134, 0.999868981, 0.999987004,
                0.981686004, 0.496752113]
        self.params = (.7928, 1.0486, -0.1336, -40.8818)
        self.bse = (.1292, .3526, .1498, 11.3979)
        self.scale = 3.3062
        self.df_model = 3
        self.df_resid = 17

    def conf_int(self):
        return [(0.5395,1.0460),(.3575,1.7397),
                (-.4271,.1599),(-63.2213,-18.5423)]
